Tuning the Kosterlitz-Thouless transition to zero temperature 
in Anisotropic Boson Systems 
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We study the two-dimensional Bose-Hubbard model with anisotropic hopping. Focusing on the 
effects of anisotropy on the superfluid properties such like the helicity modulus and the normal-to- 
superfluid (Berezinskii-Kosterlitz-Thouless, BKT) transition temperature, two different approaches 
are compared: Large-scale Quantum Monte Carlo simulations and the self-consistent harmonic 
approximation (SCHA). For the latter, two different formulations are considered, one applying 
near the isotropic limit and the other applying in the extremely anisotropic limit. Thus we find 
that the SCHA provides a reasonable description of superfluid properties of this system provided the 
appropriate type of formulation is employed. The accuracy of the SCHA in the extremely anisotropic 
limit, where the BKT transition temperature is tuned to zero (i.e. into a Quantum critical point) 
and therefore quantum fluctuations play a dominant role, is particularly striking. 
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I. INTRODUCTION 

In recent years, much progress has been made in ultra- 
cold atoms loaded in optical lattices^— Several exper- 
imental groups have demonstrated the large tunability 
of such systems by driving them from a superfluid to a 
Mott insulator phase (and vice versa) in various dimen- 
sions and lattice geometries By varying the laser 
intensity along one or several directions, experimental- 
ists can control the hopping anisotropy for the atoms 
in the optical lattice 4i§i2iifl Thus, some of these exper- 
iments have started to explore the fascinating behav- 
ior of ultracold atoms confined to low dimensions 4i5i2iifl 
This control makes it also possible to study dimen- 
sional crossover s 4 ' 1 as well as a wide range of other 
phenomena; 4 -^ which are also relevant for the under- 
standing of complex solid-state systems such as layered 
superconductors 1 ^ - — and anisotropic magnetic materi- 
als^S 2 - 

Indeed, matter in low dimensions is known to dis- 
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FIG. 1: Schematic Phase Diagram for the Bose-Hubbad 
model with hopping anisotropy ratio in two dimensions. 



play a wide range exotic properties, which are other- 
wise hard to come across in three dimensional systems. 
These include fractionalization of quantum numbers, 23 
critical states lacking long range order ^ 8 ' 10 ' 24 and topo- 
logical phase transitions that cannot be characterized by 
an order parameter such as the Berezinskii-Kosterlitz- 
Thouless (BKT) transitionJ^^ The question of how 
these exotic properties evolve as low-dimensional sys- 
tems are coupled and become, by virtue of the cou- 
pling, higher dimensional systems has attracted a great 
deal of experimental and theoretical attention in recent 
years 4- 11 - 14 - 21 - 22 - 24 - 27 - 29 

In bosonic systems, such as ultracold gases of bosonic 
atoms or molecules, as well as in anisotropic magnetic 
materials^ 4 , a theoretical analysis of the dimensional 
crossovers and other interesting phenomena such as de- 
confinement transition o 13 ' can be carried out through a 
combination of perturbative renormalization-group (RG) 
and mean-field theory (MFT) approaches. MFT assumes 
the existence of a Bose-Einstein condensate and it is 
expected to be a reliable description of the anisotropic 
superfluid phase only if the crossover takes place from 
one to three dimensions. Howecer, when trying to de- 
scribe the crossover from one (ID) to two dimensions 
(2D) at finite temperatures, MFT breaks down because 
bosons in two dimensions fail to condense at all temper- 
atures except at T = 0. Nevertheless, under such con- 
ditions a qualitative understanding of the properties of 
the anisotropic superfluid phase can be still obtained by 
means of perturbative RG and variational methods, as 
we shall demonstrate below. However, an independent 
check of these approximated methods is still required. 

It is worth noting that the superfluid properties ap- 
pear to be strongly dependent on the system dimension- 
ality—" 3 ^— Experimentally, a superfluid response has 
been observed at finite temperatures in both two— and 
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one-dimensional 30 interacting boson systems, which lack 
of a Bose-Einstein condensate. However, the origin of 
the superfluidity in these two cases is very different : 31 ' 34 
Whereas in 2D the superfluid response is essentially a 
thermodynamic phenomenon that is quantified by the 
helicity modulus, 31 in ID it is a dynamic property as the 
helicity modulus vanishes at all temperatures (the helic- 
ity modulus at zero temperature is obtained by taking the 
T — >• limit after taking the thermodynamic limi t 32 ' 34 ). 
The vanishing helicity modulus of the ID Bose fluid is 
in stark contrast with the universal jump exhibited by 
the helicity modulus across the BKT transition. As it 
will be discussed below, making the hopping amplitude 
in one direction vanishingly small, drives the BKT tran- 
sition temperature to the absolute zero (at T = 0) and 
the transition thus becomes a 3DXY quantum critical 
point (QCP) at the end of a line of classical 2DXY crit- 
ical points (cf. Fig [lj . This critical line separates the 
anisotropic normal and superfluid phases. Therefore, it 
can be argued that the helicity modulus vanishes in ID 
Bose fluids because these fluids share the same superfluid 
properties as the normal fluid phase in the limit of van- 
ishing anisotropy ratio. This can be seen by noticing the 
the helicity modulus of the ID Bose fluid can be obtained 
by approaching the ID limit along a finite temperature 
(i.e. T > 0) trajectory (cf. Fig.HJ, which means no criti- 
cal line is crossed for sufficiently small starting anisotropy 
ratio. On the other hand, approaching the ID limit along 
the T = line necessarily implies reaching the QCP first, 
which is a thermodynamic singularity (cf. Fig [lj . 

Indeed, the variety of phenomena that can be studied 
in anisotropic bosonic systems is very wide . 13 ' 14 ' 24 In this 
work, we focus on understanding the properties of the 
anisotropic superfluid phase that can be realized in e.g. 
two-dimensional optical lattices with hopping anisotropy. 
However, our results can be also of relevance to much 
more complex solid-state systems, such like anisotropic 
magnetic insulators . 21 i 22 i 24 In particular, we are inter- 
ested in understanding how the hopping anisotropy af- 
fects the properties of the superfluid phase (i.e. the he- 
licity modulus) and the Berezinskii-Kosterlitz-Thouless 
(BKT) transition temperature from the superfluid to the 
normal fluid phase. As the BKT transition temperature 
is tuned towards T = by the hopping anisotropy, the 
importance of quantum fluctuations is enhanced. Thus, 
we expect that this will lead to important renormal- 
ization effects on the parameters of the effective 2DXY 
model that describes the line of classical critical points. 
Below we shall rely on the self-consistent harmonic ap- 
proximation (SCHA) to various limits of the quantum ro- 
tor model to estimate such renormalization effects. The 
results of the calculations based on the SCHA for the crit- 
ical temperature and the helicity modulus will compared 
with Quantum Monte Carlo simulations. 

Of course, the effect of quantum fluctuations is en- 
hanced not only by the anisotropy but also by the inter- 
particle interaction, which can drive a quantum phase 
transition from the superfluid phase to a Mott insu- 
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FIG. 2: The complete phase diagram of an anisotropic Hub- 
bard model as a function of anisotropy ratio, t y /t x , interaction 
strength, pot a /U, and temperature, T/U. Note that we keep 
U to be a constant and vary other quantities for convenience. 
This is calculated by SCHA as described in the text. The red 
line guides zero temperature transitions (from superfluid to 
Mott insulator phase), while the black line is calculated by 
QMC at very low temperature (T/U — 0.01) for comparison. 
Here we set chemical potential n/U = 0.375(the mean boson 
occupation pa ~ 1) to get the QMC results. 



lator phase at integer fillings. In various dimensions, 
such superfluid-to-Mott insulator transition has been ex- 
tensively studied both experimentally and theoretically, 
mainly in isotropic systems 1 ^—. However, the combined 
effect of anisotropy and inter-particle interactions in en- 
hancing the quantum fluctuations and destroying super- 
fluidity in two dimensional Bose systems has not been 
studied so far. In this paper, we shall show how both 
quantum and thermal fluctuations can be treated on 
equal footing in the study of anisotropic superfluid in 
a 2D optical lattice. Our results can be summarized in 
the phase diagram shown in Fig. [5J 

The outline of this article is as follows: In section [H] 
we introduce the relevant low energy models that we use 
to describe the anisotropic Bose-Hubbard model (BHM) 
in 2D. Several analytic and numerical approximations to 
the BHM are also discussed there. We also discuss the 
problem of how to estimate the "phase-stiffness" param- 
eters of the 2DXY model that describes the critical line 
of BTK transitions separating the normal and superfluid 
phases at finite temperature (cf. Figs. lll2[ ). In section Hill 
the effects of thermal fluctuations and interaction on he- 
licity modulus from small to intermediate anisotropy are 
discussed. In section IIVI we explore these renormaliza- 
tion effects in the extremely anisotropic regime, as well 
as the behavior of the BKT critical temperature. The 
conclusions of this work can be found in section [V] The 
appendixes contain the technical details of SCHA calcu- 
lations in the various limits of the quantum rotor model. 



II. MODELS AND METHODS 

The anisotropy of the single particle tunneling can be 
easily realized in an optical lattice by using different laser 
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intensities for the standing waves in the x and y direc- 
tions. In the limit of a deep lattice, in which essentially all 
particles reside in the lowest Bloch band, the system can 
be described by the single-band Bose Hubbard model: 
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where 6 J is bosonic creation operator on site Rj, and 
m = b\bi is boson ocupation operator. The tunneling 



amplitude is t{ 



tyj, if Rj 



R i 



(Rj — Rj = ±ay) with a being the lattice constant. U 
and n are the on-site interaction and chemical potential, 
respectively. 



A. The XY Model and the Self-consistent 
Harmonic Approximation (SCHA) 

Deep into the superfluid phase, the low-temperature 
behavior of the system is largely dominated by phase 
fluctuations. For a sufficiently large values of U and for 
the bare anisotropy ratio parameter 770 = Jy/Jx ~ 1, we 
can represent the boson operator as 6j = V Po + $Pi e lSi > 
where po is the mean boson occupation, while 9i and 
5pi(<^. po) describe the phase and density fluctuations 
at site Hi, respectively. After integrating out the den- 
sity fluctuations, the partition function (Z) of the sys- 
tem can be written as a (Feynman) functional integral, 
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is the two-dimensional 0(2) quantum rotor model, with 
jfj = potij being the bare Josephson coupling and 
P = 1 /T the inverse of the absolute temperature in units 
where the Boltzmann constant fc# = 1. 

At sufficiently high temperatures, the imaginary time 
(r) dependence of the phase 9i can be neglected and the 
model in Eq. @ becomes the classical ferromagnetic XY 
model: 



{id) 



jfj COS 



(6i - Oj). 



(3) 



In 2D this model has two distinct phases: in the high tem- 
perature regime, the orientation of the rotors described 
by the phase 9i is disordered. The phase correlations 
are short ranged, i.e. gij = {e iei e~ l9i ) ~ e~\ Ri ~ Rj \^( T \ 
where the correlation length £(T) a. Such behavior 
corresponds to a normal phase. On the other hand, in 
the low temperature regime, the phase correlations de- 



cay algebraically, i.e. 



IRj - R 7 -|- Q(T) , where the 



exponent a(T) is finite and related to the thermody- 
namic phase stiffness. This behavior corresponds to a 
superfluid phase exhibiting quasi-long range order. The 



latter implies the absence of a Bose-Einstein condensate, 
but the finite phase stiffness (J°) means that the sys- 
tem can sustain superflows at all temperatures below the 
Berezinskii-Kosterlitz-Thouless (BKT) temperature, T c . 
Above such temperature, vortices and anti-vortices un- 
bind and destroy the superfluid properties of the system. 
The vortices (anti-vortices) are singular configurations 
of the phase 9i, where the latter winds out by positive 
(negative) integer multiples of 2tt around a discrete set 
of points on the plane. 

The picture described above relies on the classical 
(high temperature) limit of the quantum rotor model, 
where the first term in Eq. @ (oc (d T 9i) )) is neglected. 
In other words, if we expand 



(4) 



where u) n = ^-n (n being an integer), the high tem- 
perature limit only takes into account the fluctuations 
of the 9i{uj n ) field for uj n = 0. However, the model in 
Eq. © is quantum mechanical, and the quantum fluctu- 
ations are described by the finite Matsubara frequency 
(i.e. uj n ^ 0) components of 9i{uj n ). The latter and 
the classical (i.e. thermal) configurations described by 
6i(oj n = 0) arc coupled non-linearly through the Joseph- 
son coupling term oc J?- cos (9i(r) — 9j(r)). At low tem- 
peratures, both quantum and classical fluctuations must 
be taken into account. This means that we must obtain 
the effective classical limit of the quantum rotor model 
by integrating out the quantum fluctuations described by 
the 9i{uj n ^ 0) components of the phase. This is espe- 
cially important for the anisotropic XY model because, 
as we drive the system towards the extremely anisotropic 
limit where t v /t x <C 1 , the BKT transition temperature 
T c tends to zero (cf. Figs. Q] and [2]). 

In order to carry out the integration of the quantum 
fluctuations, we shall rely upon the self-consistent har- 
monic approximation (SCHA)i 13 ' 35 Thus, we shall as- 
sume that, below T c , the quantum rotor model of Eq. [2] 
can be approximated by an anisotropic Gaussian model: 
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dr 



E 



(d T t 



2U 



E <h 



<hj> 



(5) 



where is the effective Josephson coupling renormal- 
ized by the interactions and the thermal fluctuations. 
The derivation of a self-consistent equation for Jy is 
given in Appendix lAl 



B. Josephson Coupled Tomonaga-Luttinger 
Liquids and SCHA 

For small values of the anisotropy ratio(i.e. for t y /t x — > 
0), it is convenient to consider a different limit of the the 
anisotropic Bose-Hubbard model introduced in Eq. (fl}. 



4 



Indeed, for t y = 0, Eq. ([lj reduces to an array of un- 
coupled ID Bose gases. For temperatures T -C t x , 
an interacting ID Bose gas is known to behave as a 
Tomonaga-Luttinger liquid (TLL)iSi Upon restoring a 
small t y (<^i t x ) coupling between the TLLs, the result- 
ing system is an array of weakly coupled TLLs, which is 
described by the following effective Hamiltonian i 14 i 24 



Hi 



CTLL 



-y 



dx 



K {dJif + K- 1 (d^if 



p 

^4 E / dx cos [0 t - 6 H 



(6) 



where v is the sound velocity, K is the Luttinger pa- 
rameter characterizing the decay of correlations, a ~ a 
is short-range cutoff, and g°j ~ 2nt y p a 2 l /v. The fields 
^d x <j)i(x) and 8i(x) describe the (long wavelength) den- 
sity and phase fluctuations of the ID interacting Bose gas 
at site i = I, . . . , L y of the array. 

In order to obtain a phase-only description, we inte- 
grate out the density fields <pi(x) in Eq.© and thus ob- 
tain the following action for the array of weakly coupled 
TLLs: 



>CTLL\Pi 



K 
2^ 



E 

i=l 





dr / dx 
n Jo 



(d T t 



/ dr 



l Jo 



dx cosf 



), (7) 



It is now possible to apply the SCHA to this model 
by approximating the non-linear Josephson coupling in 
Stll[0i] by a Gaussian coupling: 



k v r 13 



dx 



gjv 



E 



dr 



{d T 0i) 



dx \ 



v{d x 6i) 



?i+l) 2 - 



(8) 



where gj is the effective SCHA coupling. It can be com- 
puted by solving the equation in Appendix [Bj 



C. BKT Transition and The sine-Gordon Model 

The advantage of the Gaussian models (either ([5]) or 
©), obtained after the application of the SCHA approx- 
imation, is that they allow for readily integrating out the 
" quantum components" of the phase field (i.e. the u n ^ 
components of 9) . We can thus obtain, in the continuum 
limit where the variation of the phase is slow over the 
scale of the lattice parameter, a classical Gaussian model 
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dxdy K x (d x 6) + K y (d y t 



(9) 



where the expressions for stiffnesses K x and K y de- 
pend on the starting Gaussian model: K x = j3J x and 
K y = (3J y , for Eq. ([5]), and K x = f3Kv/(a-n) and 
K y = Pgj/(cnr), for Eq. Interestingly enough, the 
role of the anisotropy in the continuum limit description 
based on ([9]) seems to be rather minor. This can be seen 
by rescaling the coordinates x — » r^l^x and y — > 77 X ^ 2 S/, 
where r\ — \J K x /K y , yielding the following isotropic 
Gaussian model: 



Sc-g[0] 



2 



dr (yey 



(10) 



where Kp = U K x K y . Note that, in a finite system, 
the rescaling also affects the system dimensions: L x — > 
LxT) 1 / 2 and L y —> LyV^ 1 ! 2 . This observation will be 
important below. 



Eq. (|10[) can be regarded as the naive continuum limit 
of Eq.(j3]) and it can only describe the (thermal) phase 
fluctuations within the superfluid phase of Eq. (fT]). Thus, 
this model can only capture the algebraically decaying 
phase correlations characterizing the superfluid phase of 
the XY model (cf. Sec. Ill A] ). However, it is completely 
unable to capture the vortex and anti-vortex unwinding 
that ultimately drives the BKT transition. 

In order to capture the possibility of topological exci- 
tations that ultimately lead to the BKT transition, we 
need to take a step back to the original XY model, ei- 
ther Eq. ([2]) or Eq. ([7]), and acknowledge that by relying 
on the SCHA, since we have neglected the possibility of 
topological configurations of the phase where the latter 
jumps by multiples of 27r from a given lattice site to a 
neighboring site. Thus, the right way to proceed would 
have been to start from the quantum XY model (or bet- 
ter, from the Bose- Hubbard model of Eq. (fT])) and, after 
integrating out the quantum components of the phase 
(and density) fields, to arrive at an effective classical XY 
model like Eq. (|3]), with properly renormalized parame- 
ters. The latter, via a duality transformation ) 36 ! 37 can 
be mapped onto the sine-Gordon model, 



S sG = I dv 



[V</>(r)] 2 2 ff (°) 



2K 



(0) 



cos 2nd>(r) 



(11) 



where </>(r) is a field that is dua l 36 ! 37 to 9(r) and g° oc 
e -E c /k B T j g j. Qe so _ ca ii ec { vortex fugacity with E c being 
the vortex core energy. The classical 2DXY and the sine- 
Gordon models belong to the same universality class, 
which means that, near the BKT transition they provide 
an equally accurate description of the long-wave length 
phenomena. For the 2DXY universality class, Nelson and 
Kosterlitz have shown^i using the renormalization group 
(RG) that, at the critical temperature for the BKT tran- 
sition, T c , the renormalized phase-stiffness (K^ ) ex- 
hibits a universal jump: 



Kf\T^T-) = - 



KfXT^T^ 



(12) 
(13) 
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The renormalized stiffness K satisfies a set of differen- 
tial RG equations, which describe the 'flow' of the sine- 
Gordon parameters (which correspond to ifg and 
at the scale of the lattice parameter a) as the system clas- 
sical degrees of freedom are coarse-grained in the vicin- 
ity of the BKT transition. Thus, RG equations deter- 
mine the long wavelength properties of the system, or, 
in other words, the phase of system: For Kp > 2/ir 
(i.e for T < T c ), the coupling of the non-linear term 
(oc cos 2TT(f>) in Eq. (fTTj) , which is responsible for the cre- 
ation of vortex-anti- vortex pairs, is renormalized down to 
zero, leading us back to the Gaussian model (cf. Eq. fTUj) 
that describes the superfluid phase, but with a renormal- 
ized value of the stiffness equal to Ki R \ On the other 

hand, when Kg < 2/tt (for T > T c ), the vortex-anti- 
vortex pairs unbind, which means that the coefficient of 
the cos27r</> term grows as the system is coarse grained. 
The unbinding disorders the system thus destroying the 
superfluidity (i.e. Kg — > 0), and thus the system be- 
comes a normal Bose fluid. 

However, it must be pointed out that the derivation of 
the sine- Gordon model from the original Bose-Hubbard 
model (cf. Eq. [T|) or the quantum XY model, Eq. @ 
is very hard to carry out in practice. The reason is 
that the integration of the w„ ^ components of the 
phase cannot be performed exactly due to the non-linear 
nature of the Josephson coupling. Thus, in this work 
we have chosen an alternative route, which involves us- 
ing the SCHA to obtain the Gaussian model with ef- 
fective parameters, K x and K y , from which we can ob- 
tain an approximation to the renormalized stiffness at 
T c : Kf{T c ) » K B [T C ) = ^ K x {T c )K y {T c ). As we shall 
show below by explicit comparison with QMC results, the 
SCHA provides a reasonably accurate estimate of the su- 
perfluid parameters even in an anisotropic Bose system 
where T c is driven to zero. Within this framework, an 
approximation to the critical temperature is determined 
from the condition that 



(14) 



Note that, since Kg is not the actual renormalized stiff- 
ness, it does not necessarily vanish for T > T c . However, 
in accordance with (|13[) we impose this fact by hand. 



D. QMC simulation on Bose-Hubbard model 

In order to validate the previously described approxi- 
mations, we have carried ab initio QMC simulations of 
the Bose-Hubbard model, Eq.(fT]), using the worm al- 
gorithm^ Earlier wor k 39 ' 40 on isotropic 2D interacting 
Bose systems has shown that this algorithm can be used 
to study the KT transition. However, as pointed out by 
Prokof 'ev and Svistunov in Ref. [ 4l| , the helicity mod- 
ulus depends strongly on the aspect ratio of the lattice 



employed in the QMC simulation, i.e. when the thermo- 
dynamic limit (L x>y — > oo) is taken by keeping L x /L y 
fixed in isotropic systems where t x — t y . As a result, the 
definition of superfluidity and its transition temperature 
can be different for different aspect ratios. The reason is 
that, as L x /L y is varied away from unity, the criticality 
of the system also undergoes a crossover from a classical 
2D XY to ID XY universality class. In the latter case, T c 
and the helicity modulus vanish. The crossover would be 
complete when we able to conduct simulations up to the 
thermodynamic limit. However, in finite systems finite- 
size effects prevent the system from completely reaching 
the ID XY fixed point. 

In this work, we focus on the effect of the hopping 
anisotropy, 770 = t y /t y =/= 1, on the superfluid properties. 
Therefore, we must first determine a physically sensible 
prescription to obtain the helicity modulus and hence 
the BKT transition temperature. To this end, in our 
QMC simulations we have chosen a value of the system 
aspect ratio, L x /L y , such that the excitation energy of 
a unit quantized flux is the same in both directions in 

the noninteracting limit, i.e. t x (yj^j — ty (ir) 01 

L x /L y — \Jt x jt y . For example, for t y /t x — 0.1, we 
use L x = 100 and L y — 32 so that L x /L y = 3.125 ~ 



Wtx/ty — vTO = 3.1622. The rationale for this choice is 
explained in what follows. 

The helicity modulus can be defined asi 42 ' 43 

2AF(tt> XtV ) 



lx,y 



(15) 



Q((fix,y / L x ^ y )^ 

where Q = L x L y is the system area and AF(<f> Xty ) is the 
free energy change due to an infinitesimal phase twist 
4>x.y applied at boundaries of the system. However, in 
a QMC simulation, the helicity modulus can be also ob- 
tained from the winding number fluctuations (W xy ) i 41 ^ 43 

T^{Wl v ), (16) 



y^x 



where (W x ) ((W y )) are the winding-number fluctuations 
along x (y) direction. 

In continuum systems, the helicity modulus, 7, can be 
related to a quantity with dimensions of density, namely 
the superfluid density p s , by means of the equation: 
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(17) 



where m is the particle mass. In lattice systems, a natural 
generalization of (|17p is obtained by replacing m by the 
effective mass, m* y ,which may be direction dependent. 

Indeed, for free particles, ^p— ~ 2t x>v , which implies 

that 



"fx.y — 2t x ^ y Ps 



(18) 



Therefore, the choice of aspect ratio L x /L y = \Jt x jt y 
means (cf. Eq. [16]) that (W 2 ) = (W^) and thus the su- 
perfluid density p s alone determines the helicity modulus 
in both directions. 
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FIG. 3: (a) shows the helicity modulus in the y direction, 7 W 
obtained by QMC as a function of the longest side, L x . We 
use rj = J°/J° = 0.02 , t x /U = 0.25 , and p ~ 0.8. From 
top to bottom lines are for T/U = 0.01 to 0.09 in steps of 0.02. 
(b) T c as a function of L x for 770 = J°/J° = 0.02. t x /U = 0.5, 
po — 0.7 for red triangles and t x /U = 0.25, po ^ 0.8 for blue 
dots. 



Next, let us assess the importance of finite size effects 
using the above choice for the system aspect ratio. In Fig. 
ED(a), we show the helicity modulus in the y direction, j y , 
as a function of the longest side, L x , where t x /U = 0.25, 
po ~ 0.8 and 770 = Jy/J x — 0.02. By the variation of 
T/U from 0.01 to 0.09, we find that, for L x > 100 and 
T < 0.05, "f y is almost unchanged. In Fig. G^b), we 
show T c as a function of L x for 770 = Jy/J x — 0.02. It 
can be seen that the variation of T c (see Sec. IIIIBI for 
an explanation of how T c is estimated from the QMC 
data) with L x is less than 5%. These results justify that 
neglecting finite-size effects on the helicity modulus and 
T c for the typical system sizes employed in our QMC 
simulations (L. x > 100). 



III. SMALL TO INTERMEDIATE 
ANISOTROPY 



A. Inside the Superfluid Phase 

We first discuss the results of SCHA for the XY model, 
which approximates Eq. @ by the Gaussian model of 
Eq. ([5]) with an effective quadratic coupling, Jij(T, U, no). 
The derivation of the equation for Jy is given in Ap- 
pendix [X] (cf. Eq. (|A12|) ). We note that the non-linear 
Josephson term in Eq.® couples all Matsubara frequen- 
cies, and therefore, in SCHA, the renormalized Jjj in 
Eq.© acquires a temperature dependence. 

From the continuum limit of the Gaussian model ob- 
tained from the SCHA (cf. Eq. [5]), the helicity modulus 
can be read off: j Xt y = 2J xy . Hence, we can also de- 
fine anisotropy ratio as 77 = J v /j x — Jy/Jx- For later 
purposes, it is also worth introducing the bare (i.e. un- 
renormalized) system parameters: 7° = 2J° = 2t XtV po 
(po is the mean lattice occupation and a the lattice pa- 
rameter) and the bare anisotropy ratio 770 = 7^/7° = 
Jy/Jx = t v /t x . In what follows, we compare the results 
of Jxyy obtained from SCHA and QMC within anisotropic 
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FIG. 4: The renormalized helicity modulus, 7^/7°, as a 
function of the interaction strength, U/J x ,. Note that, on 
both figures, we keep constant the values of U, T, and 
f]o — Jy/J x = 0.5 and 0.1, and change J x . Blue solid (red 
dashed) lines represent the results in the a — x(y) directions. 
For comparison, we also show the numerical results obtained 
by QMC in filled circles and triangles together. See the text 
for more details of comparison. 



superfluid (SF) phase. 

In Fig. U we show the ratio of the renormalized to 
the bare helicity moduli, 7^/7^ as a function of the 
interaction strength, U/J x . The bare anisotropy ratio 
parameter is chosen to be 770 = 0.5 ( Fig. IDJa)) and 
7/0 = 0.1 (Fig. Fig. IDJb)). Here we keep both U and 
T constant but change J x and J° in order to compari- 
son with QMC data more easily. As expected, increas- 
ing the strength of interactions, that is, increasing U /J x , 
suppresses superfluidity as that both 7^ and j y decrease. 
Note that, within the SCHA, even a weak interaction can 
have a strong effect on the renormalized helicity modu- 
lus, J x ,y Indeed, when the interaction is larger than 
a critical value, the helicity modulus drops to zero in 
both directions discontinuously, and the system becomes 
a normal fluid without phase stiffness. This is a feature of 
the SCHA, which wrongly predicts the interaction-driven 
transition between the SF and the Normal fluid (which at 
T = corresponds to the SF to Mott insulator quantum 
phase transition) to be of first order. 

In the same figure, we also show numerical results of 
our QMC simulation for comparison. We can see that, 
although the ratio of the renormalized to the bare helic- 
ity modulus obtained from QMC exhibits qualitatively 
the same behavior as the SCHA, it does not show strong 
renormalization effects predicted by the SCHA at small 
U/J x . Furthermore, at larger U/J x , both 7^/7° and 
Jy/jy vanish rather smoothly. 

In order to better understand how finite temperature 
and interactions influence the anisotropy ratio of the he- 
licity modulus, we show in Fig. [5ja) and (b) the renor- 
malized helicity ratio, 77 = 7^/7^ vs. the bare one 
Vo = lyllx = ty/tx- Results obtained both from the 
SCHA and QMC are shown together for comparison. We 
see that, at low temperatures (T/U — 0.01, Fig. 0a)), 
when the system is deep in the superfluid phases, the 
anisotropy is barely renormalized, i.e. 7/ ~ 770 and indeed 
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FIG. 5: Ratio of the renormalized helicity moduli, rj vs. its 
bare (unrenormalized) value, rjo, for (a) T/U = 0.01 and 
(b) T/U = 0.5. We choose different values of J°/U: (a) 
J°/U = 0.325 and J°/[7 = 0.188 and (b) J°/U = 0.275 and 
Jx/U = 0.163, and vary 770- Note that in (b), for 770 J$ 0.5, 
the temperature is higher than the BKT transition tempera- 
ture and therefore both 7^ and y v vanish. We also show the 
results of our QMC simulations for comparison. Both QMC 
and SCHA yield results in excellent agreement, suggesting 
that the anisotropy ratio of the helicity modus is barely renor- 
malized by interaction and finite-temperature effects. This is 
consistent with the SF phase being described by an isotropic 
Gaussian field theory, as discussed in Sec. Ill CI 



our QMC results agree well with the SCHA predictions 
for 77. Interestingly, this result holds also true at much 
higher temperatures (cf. Fig. |5Jb)), except for the fact 
that, for small values of 770 the system becomes a nor- 
mal gas (i.e. the temperature used in the simulation 
T/U = 0.5 is larger than the BKT transition temper- 
ature, T c , for these highly anisotropic systems). This is 
because as T c of an anisotropic superfluid (770 < 1) be- 
comes smaller, the renormalization of the helicity ratio 
also becomes more significant near the phase transition 
boundary. The agreement between SCHA and QMC re- 
sults is very good. 

Thus, we find that, although the SCHA and QMC yield 
different values for renormalized helicity moduli, j x and 
7 y , QMC shows that the renormalized anisotropy ratio 
77 = "fy/^fx is barely affected by interaction and/or tem- 
perature effects. This is consistent with the SF phase 
being described, in the continuum limit, by an isotropic 
Gaussian field theory (cf. Eq. [TO] in Sec. Ill C[) . which is 
also correctly captured by the SCHA. 



B. Near the BKT transition 

In Fig. [6l we show the helicity modulus (proportional 
to superfluid density cf. Eq. [17]) as a function of tem- 
perature. The bare single particle tunneling amplitude 
is t x /U = 0.5 and t y /U — 0.25 respectively, and the fill- 
ing fraction is po ~ 0.63. In Fig. [U(a), the results of 
the helicity moduli obtained from the SCHA and QMC 
are compared. We thus see that, compared to the QMC 
results, the SCHA overestimates the temperature depen- 
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FIG. 6: The helicity modulus as a function of temperature 
for a fixed interaction U. The bare tunneling are t x /U = 0.5, 
ty/U — 0.25, and the density is po ~ 0.63. (a) shows 
results including interaction renormalization within SCHA, 
compared with the QMC results in dots and in triangulars. 
(b) shows the winding number fluctuation in QMC, and Kp 
obtained analytically from the SCHA to XY model, as a func- 
tion of temperature. The horizontal dashed line indicate the 
universal number, 2/tt. The intersection of the curves (dots) 
and the horizontal lines gives the T c in SCHA (QMC), marked 
by arrows. 



dence of the helicity modulus in both directions roughly 
by a factor of one point five. 

In Fig. [D(b), we show how the BKT transition temper- 
ature T c is determined from both the analytical results 
of SCHA and the QMC data. In the case of the SCHA, 
we compute the phase stiffness as discussed in Sec. Ill CI 
i.e. from Kp = /3 •»/ J x J y , where J x and J y are solutions 
to the SCHA equations for given T and U,J X ,J X val- 
ues. Hence, T c is found by varying the temperature until 
Kp = ± (cf. Sec. ESI). 

As to the QMC data, T c is obtained as follows: As 
anticipated in Sec. Ill Dl by choosing L x /L y — \Jt x /t y 
we find that the winding number fluctuations (red dots 
and black triangles in Fig.[6jb)) in the x and y directions 
essentially coincide. Furthermore, the (W% y ) show kink 
at a temperature, which is essentially equal to the one 
obtained by requiring that T c (cf. Sec lII CP : 



K$ MC (T C ) = l^ Jx (T c ) 7y (T c ) 



J(W*(T e ))(W*(T c )) ~ (Wl l 



- (19) 

IT 



where Eqs. (|T4|) and (fT6|) have been used. In Fig. [6jb), 
we have indicated the universal value of — by a hori- 
zontal line. As explained in Sec. Ill CI in the SCHA, we 
assume that Kp vanishes for for T > T c . However, in 
the QMC calculations, finite-size effects round off the ex- 
pected thermodynamic-limit discontinuity of K ( ^ MC at 
T = T c . Yet, as discussed in Sec. Ill Dl the value of T c 
estimated from the kink in the Monte Carlo data is con- 
verged for system sizes that we used (L x > 100). Finally, 
the comparison of T c as obtained from QMC and SCHA 
is shown in Fig. [S] and will be explained in more detail 
further below. 
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FIG. 7: The Luttinger parameters K and v as a function of 
the temperature, T/U. For the ID Bose-Hubbard model, we 
fix L x — 150 for both plots, and (a) t x /U = 0.5, t y /U = 
0,p ~ 0.7 and (b) t x /U = 0.25, t y /U = 0, p ~ 0.8 respec- 
tively. 



IV. LARGE ANISOTROPY 
A. SCHA for the Coupled TLLs 

To begin with, let us note that, for the ID Bose- 
Hubbard model, the Luttinger parameters K and v that 
determine the properties of the TLLs in the decoupled 
limit (cf. Eq. [5] for g°j = 0) cannot be analytically ob- 
tained for general lattice fillings and values of U/t x (Eq.[T] 
for t y = 0). Thus, in order to extract the Luttinger liquid 
parameter, K , and sound velocity, v, we have carried out 
additional QMC calculations for the ID Bose-Hubbard 
model to extract these parameters. Using the relations 
v/K = 1/ttk and vK — irL. x T(W x ), where k = dp/ dp, is 
the compressibility and (W x ) is the winding number fluc- 
tuation along the x direction for T/U -C l- 44 In Fig. [7Jthe 
numerical Luttinger parameters K and v as a function of 
the temperature are shown, T/U, for a large large size 
of the ID system of L x = 150. The parameters charac- 
terizing one (decoupled) TLL correspond to the extrap- 
olation of this results to very low temperature. Thus, 
for T/U = 0.005 we find K ~ 2.77 and v ~ 0.77 for 
t x /U = 0.5, po ~ 0.7 in Fig. [TJJa), and K ~ 1.91, v ~ 0.53 
for t x /U = 0.25, po ~ 0.8 in Fig. [7J». 

Next, we describe the result of applying the SCHA 
to the system of coupled TLLs. Compared to the case 
of small anisotropy discussed above (Eq. IA12|) . in this 
case only J y is renormalized, and all the interaction de- 
pendence of J y enters through the Luttinger parameters. 
However, as discussed in sec. Ill CI the system of coupled 
TLLs at finite temperature also belongs to the 2DXY uni- 
versality clas o 13 ' 14 . Thus, the BKT critical temperature 
can be found from the equation: 



K B =2 



y/KvJy(T c )/2ir 



(20) 



In the SCHA calculations, we have chosen the short- 
distance cut-off such that Kvclq /2tt ~ J x , when compar- 
ing the TLL-Gaussian model of Eq. © with the Gaussian 
model in Eq. ((3J). 



0.8 

0.6 

y"0.4 
H 

0.2 



(a) 



0.0 



r 



• SCHA(XY) 
I SCHA(TLLs: 
■ QMC 



(b) 



p 



I SCHA(XY) 
i SCHA(TLLs) 
■ QMC 



0.0 0.2 0.4 0.6 O.i 
*/0 



1.0 0.0 0.2 0.4 0.6 0.8 1.0 

1o 



FIG. 8: T c as a function of the bare anisotropy ratio, 770 = 
J°/J° = t y /t x for t x /U = 0.5, po ^ 0.65 ± 0.05 (J° ~ 0.325) 
(a) and t x /U = 0.25, p ~ 0.75 ± 0.05 {J°/U ~ 0.1825) (b) 
The blue triangles are the results obtained by SCHA to the 
XY model. The black dots correspond to the results obtained 
by SCHA to an array of coupled TLLs. The Red squares are 
the QMC data. T c is determined by the methods discussed 
in Section l!!! Bl On both panels, the green dash curve is a fit 
to the scaling behavior of Tc with the bare anisotropy ratio 
770 yielding T c /U ~ 0.837»7g' 55 (a) T c /U ~ 0.448r;g- 575 (b). 



In Fig. [8] (black dots), we show the BKT critical tem- 
perature, T c computed using the SCHA, and QMC as a 
function of the bare anisotropy ratio. As discussed above, 
T c goes to zero gradually as the bare anisotropy ratio 770 
becomes larger, reflecting the fact that there is no super- 
fluid phase transition at finite temperature in ID system. 
From both panels in Fig. [SI it can be seen that the SCHA 
to XY model provides a reasonably good description of 
T c (compared to the QMC results) for 770 ~ 1 and weak 
interactions (Fig. [8ja)), but it deviates from the QMC 
results for stronger interactions (Fig. [BJb)) and small 770. 
On the other hand, the results obtained by applying the 
SCHA to an array of coupled TLLs are found to be closer 
to the QMC results for T c in the large anisotropy regime 
(i.e. small 770). These results are consistent with the 
expectation that the SCHA to the XY model should be 
more accurate in the small anisotropy regime, whereas 
applying the SCHA to an array of coupled TLLs becomes 
a better approximation in the limit of large anisotropy. 



B. RG scaling for critical temperature 

Besides of the numerical calculations of BKT critical 
temperature, from our QMC data we can also extract the 
scaling behavior of Tc with anisotropy ratio 770 ■ This can 
be compared with the results obtained by the renormal- 
ization group flow of the Josephson coupling in Eq ([7)1. 
which described by the differential equation ; 13 ' 1 



dg.j 
d£ 



= 2 



1 

2K 



gj- 



(21) 



where the flow parameter £ = ]n(a(£)/ao) = 
ln(A(0)/A(i)), with a{£) = a e e . Since K € [l,+oo) 
for the Bose-Hubbard model is far from the critical point 
K* = 1 /4, we can neglect the renormalization of K and 
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treat it as a constant ! 13 ' 14 Therefore, the solution to (f2"Tj) 
reads gj(£) — g,j(0)e^ 2 ~ trtW. To complete the solution, 
we need to recall that the bare (energy) cut-off A(0) ~ t. x , 
and gj(0) — In order to estimate of the crit- 

ical temperature at which the system will enter the SF 
phase, we note that, at finite temperatures, the RG flow 
is cut off at the scale A(£) ~ T, and gj(T = T c ) ~ 1. 
Hence, provided (|2"Tj) provides accurate description of the 
e entire flow (i.e. for small enough gj(0)), we have 



(22) 



where C is a prefactor that depends on microscopic de- 
tails of the model, and can be obtained by fitting above 
scaling law to the QMC results. It is worth noting that 14 
the same scaling law for T c /t x with r/o can be also ob- 
tained using mean-field theory, i.e. by assuming that 
(e* e "(°)) = (j>o(T). However, as discussed in the Introduc- 
tion, strictly speaking mean-field theory is inapplicable 
in two dimensions due to the lack of BEC at finite tem- 
peratures. 

In Fig. [8]we use the values of the Luttinger parameters 
obtained earlier from QMC simulations of the ID Bose- 
Hubbard model (K ~ 2.77 and v ~ 0.77 in Fig. [gfa) 
and K ~ 1.91, and v ~ 0.53 in Fig. 1Kb)) to fit the 
scaling of T c . In particular, the value of K completely 
determines the exponent of the scaling law (cf Eq l2"2l 
and thus, the only free parameter is the prefactor C 
The fit yields T c /U ~ 0.837ryg' 55 for the data on Fig 
UKa) and T c /U ~ 0.448t# 575 for the data on Fig. HJb) 
Using the TLLs parameters obtained from T/U — 0.005 
the predicted T c of SCHA to TLLs are close to QMC 
calculations at large anisotropy regimes. 



V. CONCLUSION 

In summary, two different approaches, the SCHA and 
QMC, reveal the highly nontrivial features of the helic- 
ity modulus and the BKT phase transition in the 2D 
Bose-Hubbard model with anisotropic hopping. These 
characteristic features simulated by QMC using a spe- 
cific system aspect ratio, L x /L y = y/t x /t v , is consistent 
with the rescaling of the effective sine-Gordon model. We 
show how the interaction and finite temperature effect 
influence the helicity modulus and find profound agree- 
ment of anisotropy of the helicity modulus given by the 
SCHA and QMC. As we drive the system towards the 
extremely anisotropic limit, the BKT transition temper- 
ature approaches to the absolute zero and the transition 
thus becomes a 3DXY quantum critical point at the end 
of a line of classical 2DXY critical points. In particular, 
through the RG scheme for the coupled TLLs, we obtain 
the scaling relation of Tc with anisotropy ratio. Em- 
ploying ultra-cold atoms in an controllable optical lattice 
opens avenues to identify our results of 2D anisotropic 
Bose-Hubbard model. 
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Appendix A: Self-consistent Harmonic 
Approximation 

To find the optimally quadratic approximation to the 
XY-model, we employ the self-consistent Harmonic ap- 
proximation (SCHA). In this approach, the action of XY 
or quantum rotor model, Eq. ([2]), is approximated by an 
anisotropic Gaussian model: 

Sam = jf *-{£^£+ E JiM-m 



<i,j> 



(Al) 



where w„ = 2nTn, and the single particle Green's func- 
tion is given by 

G-^wJ = ^+^8J Q sin 2 (fc Q a Q /2). (A2) 

a 

Next, we make use of Feynman's variational principle, 
which states that: 

F = ~\nZ< F[G V ] = F V + ^(S[9} - S G [9]) V , (A3) 

where () v denotes the average with respect to Sq[6] and 
S is the XY model action. Since 

- / D6e- S ^ = JJ G v {\w n )-^ , (A4) 



the first term of F[G„] is: 

F v = - go £ \nG v (k,u n ). 

The remaining contributions to .F[G„] are 

(Sx y [e]-s G [9]} v 

ri 3 



(A5) 



</ <M£^(Mi) a - E 2J°.cos(0,-0 J )}}„ 

i <i-,j> 

-{Sg[6])v 

£ -^rG^(k,o; n ) + (S cos ) v - const. (A6) 



with 



(Scos)v 



/ dr £ 2J^{cos{9 l ~e 3 )) v . 
Jo 



<h3> 
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Hence, 



Using the Matsubara sum 



(cos(6>,: — 9j)) v = Re 
= Re 



e -i((e 1 -e J ) 2 )„ 

Gv(Xi-Tjfl)-G v (0,0) 



(A7) 



by the cumulant expansion. Here 



G„(w n ,k), 



is the single particle Green's function in real space. 
Therefore we have, 



(S cos ) v = / dr V 24 Re 

= -/»E E 2J ° Re 



3 G„(i-i-rj,0)-G„(0,0) 



,G„(t,0)-G„(0,0) 



i t=S-,a„ 



-,80 E 2 J°Re [e W (e*--- l )G.(k, UB ) 

a 

(A8) 



Upon combining above results and finding the extrema 
of F[G V ], i.e. 



6F'[G V ] 
6G v Qc,u) n ) 



= 0, 



(A9) 



we find 



1 



G„(q,w„) 
2 



^G„(k,u„) 



[7 



2^ C ° th( — < 



(All) 



we conclude that 



l)^C0,h(^).(A12) 



Note that Wk = 2\/2t7y J x sin 2 (fc x a/2) + J y sin 2 (fc a a/2) 
is the phonon (Bogoliubov) excitation energy. 



^ + 8 ]T J° sin 2 (g Q a Q /2) e W -i)G„(k,u„ 
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Appendix B: SCHA for Coupled TLLs 



Applying the methods of previous section to the action 
of Eq. ([7]), the following equation for the renormalized 
parameter J y (= gjv/-Ka^) is obtained: 



^ = i^E-^^^-a-). ( B1 ) 



^E Jq sin 2 (q Q a Q /2). 



(A10) where w k = 2^v 2 (k x /2) 2 + 2 ^ sin 2 (fc^/2). 
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